East European sedimentary basins long heated by a fading mantle upwelling

A strong negative anomaly of seismic wave velocities at the core-mantle boundary (the Perm Anomaly) beneath the East European platform is attributed to the remnant of a deep mantle upwelling. The interaction between the upwelling and the East European lithosphere in the geological past and its resulting surface manifestations are still poorly understood. Using mantle plume modelling and global plate motion reconstructions, we show here that the East European lithosphere is likely to have been situated over the weakening Perm Anomaly upwelling for about 150–200 million years. As the East European platform moved above the Perm Anomaly in post-Jurassic times, the vertical tectonic movements recorded in sedimentary hydrocarbon-rich basins show either hiatus/uplift or insignificant subsidence. Analytical modelling of heat conduction through the lithosphere demonstrates that the basins have been slowly heated for a long time by the Perm Anomaly upwelling, creating suitable conditions for hydrocarbon maturation. This suggests a profound relationship between mantle plume dynamics, basin evolution, and hydrocarbon generation.

extinction event [1][2][3] .Fluid dynamic experiments suggest that such a plume likely originates from a thermal boundary layer at the core-mantle boundary 4 .
Reconstructing the Paleozoic mantle remains speculative due to the challenges in determining the deep mantle sources, which generate upwelling, and the paths of the upwelling within the convecting mantle.Establishing reliable Absolute Plate Motions (APMs) is further complicated by uncertainties related to tectonic deformations and continental arrangements, particularly, the complex kinematics of the North Atlantic bordering continents.Consequently, various models have been proposed including those linking the Siberian traps to the Iceland hot spot 5 .
Meanwhile, the paleogeographic and palinspastic reconstruction based on paleomagnetism 6 and true polar wander-corrected paleomagnetic reconstructions 7,8 have provided evidence of a possible connection between PA and the Siberian trap province, if the source of the PA upwelling remained reasonably fixed in space.According to these reconstructions, the Siberian trap eruption site was located relatively close to PA by the end of the Permian times.Considering the mobility of the PA plume, Flament et al. 9 proposed that the PA can be linked to the Emeishan large igneous province (LIP).However, Torsvik & Domeier 10 dismissed this hypothesis citing the stability of large low-shear-velocity provinces (LLSVP), the location of the Emeishan LIP eruption site near the equator 260 My ago, and its association with the JASON LLSVP (see Fig. S1).Similarly, it is plausible to reject the notion that the PA mantle plume is linked to the Devonian Kola-Dnieper LIP (KDLIP) 11 .According to the paleogeographic and paleomagnetic reconstruction 12 , the eruption sites of KDLIP were situated near the equator approximately 360 My ago and were potentially associated with the TUZO LLSVP.Therefore, KDLIP was not connected to the present PA (Fig. S1).

SD2. Magmatism of the East European platform since the Triassic times
If PA is considered to be a remnant of the old mantle plume and this plume is likely to be linked to the Siberian LIP, major magmatic events would have been associated with the Permian-Triassic boundary (the time of the Siberian traps).In the beginning of the Early Triassic, magmatism was still active in West Siberia 13 .In the Late Triassic (225-207 My) magmatism was manifested as tholeiitic and calc-alkali basalt eruptions in the Turgay basin (with the geographic centre of the basin to be around 50°N, 66°E) located to the southeast of the Urals 14,15 .However, there have been no substantial indications of volcanism recorded in the East European (EE) lithosphere since the Early Jurassic times.A reason for that is that the EE lithosphere was thick enough to allow for magmatic eruptions.The thickness of the lithosphere exceeds 150 km in the studied area: this is visible in the S-wave velocity images of the EE lithosphere (Fig. S2) plotted using the data from Shapiro & Ritzwoller 16 , and in the lithosphere-asthenosphere boundary maps based on the data from Hoggard et al. 17 (Fig. S3) and from Artemieva 18 (Fig. S4).The absence of magmatic events in the studied region of the EE platform can be explained by the thick (>150 km) lithosphere 19 .

Fig. S2.
Six cross-sections of S-wave velocities illustrate pronounced positive anomalies of the velocities associated with the cold and thick lithosphere of the East European platform.The cross-sections (presented as a percent deviation from the velocity in the 1-D model ak135) have been plotted using the data from the model CU_SRT1.0 (ref.16) and the interactive plotting tool (http://ciei.colorado.edu/~nshapiro/MODEL).Black contours outline the persistent velocity anomalies, and small black dots show the earthquake locations 22 .The white shaded areas in the cross-sections correspond to the studied region.World maps above the cross sections illustrate their location (red lines).Fig. S3.The lithosphere-asthenosphere boundary (LAB) derived from the S-wave seismic velocity model 23 .Data are taken from 0.5°×0.5°grid model by Hoggard et al. 17 .Numbers on the isolines present the LAB depth here and in Fig. S4.

Fig. S4.
The LAB defined as the intersection of the mantle geotherm with the 1300°C mantle adiabat.Data are taken from a 1°×1° grid model by Artemieva 18 .Bold lines mark the area where the LAB depth is more than 150 km.

SD3. Global plate motion models
Three models of the global lithospheric plates motion from 250 million years (My) to the present have been developed for the last decade 8,24,25 ; hereinafter we refer to the models as Se2012 (ref.24), Mü2016 (ref.25) and Ma2016 (ref.8).The models attempted to reconstruct the APM relative to the Earth's deep interiors and used and/or discussed the APMs from the present to 70 My -~120 My were earlier developed [26][27][28] .These APMs link to the global plate circuit through Africa, as the African continent has been surrounded by mid-ocean ridges for about 170 My from the present.
The Se2012 model was based on the moving Indian-Atlantic hotspot model 27 for times earlier than 100 My, while the Ma2016 and Mü2016 models used the Torsvik et al.
model 28 (reasons for that were discussed in papers 8,25 ).The Torsvik et al. model 28 was preferred due to the more reasonable trench migration and net lithospheric rotation (NLR) amount implied.Indeed, their analysis suggests that APM models, in which Africa drifts NE during the 130-70 My period, are most likely in terms of trench migration rates and NLR, and also satisfies more accurately the slab remnant constraints proposed by the APM based on subductions alone 29 .A westward shift of ∼14 deg. at ∼130 My compared to the hotspot APM models falls within the range of possible fits to seismic tomography, whilst also yielding reasonable NLR and subduction zone kinematics 25 .
Although reconstructions with respect to mantle are quite speculative for the Mesozoic and Palaeozoic times, three maps of global plate motion were proposed until 200 My 24 , 231 My 25 , and 250 My 8 .We recall that palaeomagnetic reconstructions relate the past configurations of continents to the Earth's spin axis without palaeolongitudes.Moreover, plate reconstructions with respect to mantle must include an estimate of true polar wander (TPW).Steinberger & Torsvik 30 argued that the TPW direction for the Cenozoic and Mesozoic times was largely controlled by the mass of the two antipodal large low shear-wave velocity provinces and assumed that these deep mantle bodies have been stable over a longer time scale.Hence, they proposed an estimate of TPW from 100 to 250 My.
To summarize, Seton et al. 24 used a combined reference frame: the O'Neill et al.
model 27 for the last 100 My and the palaeomagnetically-derived TPW corrected reference frame 30 from 100 My to 200 My.Müller et al. 25 discussed extensively the present to 130 My APMs as mentioned above and finally used the Torsvik et al. model 28 .A close inspection of their table for the part of 130 My to 230 My is nevertheless necessary to understand what they have done: they used a combination of a rotation due to two interpolated longitudinal shift-TPW corrections, and an additional correction using the APM rotation poles 30 .Finally, Matthews et al. 8 also used the Torsvik et al. model 28 and computed an APM between 100 My and 250 My by using the sum between three Euler poles derived from the palaeomagnetic poles 31 , using TPW correction 31 and the same pole rotations as in Müller et al. 25 .These models can be inspected using the GPlates (https://www.gplates.org).Also, GPlates gives the Euler rotation parameters describing both the relative and the absolute motion of lithospheric plates with respect to what is thought to represent a mantle reference frame.We note that the plate circuit between Africa and Laurasia is identical in the two models 8,25 , while it offers only minor differences with the one used by Seton et al. (ref.24).To derive tracks of the PA upwelling related to the Se2012, Mü2016, and Ma2016 models, we have used the set of parameters presented in Table S1, where λ is the latitude, and φ is the longitude.The reconstruction of the PA's position relative to Siberia approximately 260 My ago, using the APM model 8 , is depicted in Fig. S5.According to this model, the PA was situated to the south of the West Siberian Basin and is likely responsible for the Permo-Triassic volcanic activity documented in drilling operations 2 .The Late Triassic basalts in the Turgay basin could potentially mark the concluding visible volcanic activity on the surface (see also SD2), in accordance with the reconstruction presented in Fig. S5.stretching of the lithosphere and over time scales much longer than predicted by the model of thermal subsidence 32 .Lobkovsky et al. 33 suggested that an intraplate compression and stressinduced lithospheric deflection 34 as well as a thermal mantle upwelling followed by eclogitization-induced subsidence 35 can play an important role in the formation of EE basins.
The tectonic subsidence analysis was performed using borehole data for the Timan-Pechora basin 36 , the Moscow basin 37 , the Pricaspian basin 38 , and the Pre-Uralian foredeep.
The location of the boreholes is presented in Fig. 3 of the main text as well as in Figs.S7, S9, and S11.To determine a tectonic subsidence, the backstripping method is employed to separate isostatic effects of sediments and water loading from those of tectonic subsidence 39 .
The backstripping method is based on assessments of the (Airy) isostatic response to sediment load and does not consider the flexural strength of the lithosphere, which means that any load -positive, in the case of deposition, or negative, in the case of erosion -influences the subsidence or uplift history, not just vertically below the load, but also laterally around the load.For a given sedimentary load, the magnitude of vertical deflection under the loads and the lateral wavelength of the deflection are controlled by the effective elastic thickness of the lithosphere, Te.To evaluate the response of the elastic lithosphere to the elevation of the topography h due to sediment load, we consider the amplitude of the deflection w of the elastic lithosphere 40 : My -1 in Pliocene (Fig. 3 in ref. 32).An increase of the subsidence rate (to about 10 m My -1 ) can be observed in some borehole data in the Quaternary times, and the mechanism for this Quaternary subsidence acceleration in the EE basins is still poorly understood.
The Pricaspian basin is situated on the south-eastern portion of the East-European (EE) platform at the northern end of the Caspian Sea.The basin is approximately 600 km across from west to east; the thickness of the sedimentary cover is more than 20 km in the basin centre 43 .The basin infill is divided into three major sedimentary sequences: subsalt strata, salt, and overburden of the salt.The subsalt sequence contains Riphean through Lower Permian strata punctuated by unconformities, and it has a complex depositional history dominated by carbonate reefs and clastic fans 43 .The salt sequence consists of Kungurian (ca. 260-258 My) salt overlain by Kazanian (ca.258-252 My) salt, which reaches a thickness of 4.5 km in the centre of the basin.The overburden of salt consists predominantly of terrigenous Upper Permian through Neogene strata 43 .The tectonic subsidence analysis 38 was performed for the central portion of the Pricaspian basin using the synthetic stratigraphic column (Fig. S6; see also Fig. 4 of the main text).The rapid subsidence in the Pricaspian basin began in the Middle Devonian, which could be associated with the Kola-Dnieper mantle plume activity since 380 My ago 11,12 .
Located to the west of the Pricaspian basin, the Dnieper-Donets aulacogen underwent stretching due to Devonian rifting 44 , basaltic volcanism in the Late Devonian times 11,45 , and an episode of rapid subsidence in the Late Devonian 32,46 .The Jurassic (208-170 My) uplift in the entire Pricaspian basin resulted in about 38 My hiatus.The tectonic subsidence analysis 37 was performed based on the 18 boreholes drilled in the central and southern parts of the Moscow Basin (Fig. S7; also see Fig. 3 of the main text).Vertical tectonic motions of the basement of the Timan-Pechora region were studied based on an analysis of 37 exploration wells 36 (Fig. S9; also see Fig. 3 S11; see also Fig. 3 of the main text) were used to analyse the regional vertical motions.The tectonic subsidence curves (Fig. S12) show no or little subsidence for more than 100 Myr.

SD5. Mantle plume diffusion modelling: Limitations, uncertainties, and sources of errors
The mathematical model of mantle plume dynamics used in this study is simple (see Methods in the main text), and many complications are omitted 47 .An increase of the viscosity from the upper to the lower mantle is not considered.The adiabatic heating and cooling as well as mantle phase transformations are not included in the model, while the phase changes retard or accelerate the ascent of mantle plumes 48 .Radiative thermal conductivity in the lower mantle may lead to even faster diffusion of plume tails 49 .
Numerical models of mantle plume evolution have several sources of errors and model To assess a sensitivity of numerical results with respect to the Rayleigh number Ra, several numerical experiments of mantle plume evolution were performed at Ra ranging from about 10 4 to 10 6 (ref.47).The experiments showed a similar pattern of the plume's decay, although the morphology of the plumes and the time of the decay were different.To reduce uncertainties in numerical models of mantle plume evolution, which are related to the initial conditions in the geological past, the present data can be assimilated to reconstruct the evolution of a present diffused mantle plume to its prominent state in the past 47,51 .With the development of modern data assimilation techniques 52  Also, in this paper, we do not analyse a detailed interaction of the continental lithosphere with the underlying mantle.It has been shown that a continental rheological structure and intra-plate stresses may influence the topography of the lithosphere above the mantle plume 54 and hence the thermal evolution of the lithosphere.

Fig. S1 .
Fig. S1.Reconstruction of large igneous provinces using a hybrid reference frame 20 and draped on the tomographic model 21 .The plume generation zones marked by thick red lines in this model corresponds to the 0.9% slow S-wave seismic wave velocity contour, and the zero contours are shown as thinner black lines.LIPs with red and green symbols reconstructed with different reference frames.SIB: the 252 My old Siberian trap province; ELIP: the 260 My old Emeishan LIP; and KDLIP: the ~360 My old Kola-Dnieper LIP.KDLIP is twocomponent event -the Kola Alkaline Province (the left violet circle), and the Dnieper-Donets aulacogen (the right violet circle), which are located on far sides of Baltica paleocontinent 11 .The position of the KDLIP corresponds to the paleogeographic/ paleomagnetic reconstruction for 360 My (ref.12).Modified after Torsvik & Domeier 10 .

Fig. S5 .
Fig. S5.Reconstruction of the PA's position relative to Siberia at the time of 260 My (Late Permian) and 210 My (Late Triassic) using the APM model 8 .Red, orange and yellow lines are the contours of the PA according to different Vs velocity anomalies.Green line marks the overall contour of the Siberian LIP.Black triangles are the oil wells in the East European platform (EEP).Blue crosses are main (present) hotspots.The green cross shows the location of Norilsk, the petrographic centre of the Siberian traps (ST); and blue crosses the present main hotspots.The absolute position of Siberia is from the APM model 8 .

ρ
is the density of the mantle; c ρ is the density of the crust; E is Young's modulus; ν is Poisson's ratio; g is the acceleration due to gravity; and λ is the topographic elevation's wavelength.The degree of the compensation max / C w w = of the topographic load can be obtained as the ratio of the deflection of the elastic lithosphere to its maximum deflection in the case of hydrostatic equilibrium.Considering the effective elastic thickness of the East European lithosphere to be Te=100 km (ref.41), E = 70 GPa, ν = 0.25, m ρ = 3300 kg m -3 , c ρ = 2800 kg m -3 , and g = 9.8 m s -2 , we can find that the topography is about 30% compensated if λ ~1000 km.The tectonic subsidence of the EE basins since the Devonian indicates that the platform underwent two major episodes of rapid downward motion: the first episode occurred within all basins in the Late Devonian, and the second episode took place within the Timan-Pechora, Pricaspian, and Moscow basins in the Late Carboniferous/Early Permian 32 .Based on the global vertical motion of the EE platform 42 , the subsidence rate of EE basins ranges from about 18 m My -1 in the Late Devonian (~360 My) to about 6 m My -1 in the Early Jurassic (~200 My).The subsidence slows down from about 4 m My -1 in the Middle Jurassic to 1 m

Fig. S7 .
Fig. S7.Locations of the boreholes in the Moscow basin Basin commenced in Precambrian time (Fig.S8).At least two episodes of rapid apparent subsidence since the Cambrian occurred within the Moscow Basin in the Devonian and in the Late Carboniferous/Early Permian, respectively.The Moscow Basin became an area of erosion since at least early Cretaceous.The Timan-Pechora Basin is situated at the north and northeast outskirts of the EE platform and is adjacent to the Urals on the east.On the north the Timan-Pechora Basin extends beneath the Barents Sea.The main structural elements of the basin include the Timan ridge and the Pechora syneclise.The Timan ridge is a major linear basement elevation along the western edge of the basin.The Pechora syneclise is identified in the sediments and occupies most of the basin area.The maximum thickness of sediments is 7-8 km in this part of the basin.It was formed above a zone of Palaeozoic subsidence at the eastern edge of Timan-Pechora during the orogenic phase in the Urals history.
of the main text) and three lithologic-stratigraphic cross-sections along the regional seismic profiles.The tectonic subsidence curves for the wells are shown in Fig.S10.These curves display a stepwise pattern involving periods of rapid subsidence followed by phases of slower subsidence.The Timan-Pechora basin underwent tectonic uplift followed by basaltic volcanism in the Middle-Late Devonian, which could be associated with the Kola-Dnieper mantle plume.The basin experienced rapid subsidence during Late Devonian time (an episode of rifting) followed by slow subsidence of the basement (an episode of thermal relaxation).In the Early Permian the basin again experienced relatively rapid subsidence diminishing in the Middle to Late Jurassic compared to the Devonian and Permian-Triassic times.Nearly all Timan-Pechora suites became an area of erosion in the Late Cretaceous and Palaeogene.

Fig. S9 .
Fig. S9.The location of boreholes in the Timan-Pechora basin

uncertainties 47 :
errors associated with the idealization of mantle plume dynamics by a set of conservation equations; errors related to the difference between the exact solution to the conservation equations and the exact solution to the discretised equations; and errors arising from iterative solutions to the model equations.In addition, we should account for errors related to unknown mantle temperature in the geological past and the uncertainties in initial and boundary conditions of the model.Particularly, the temperature at the core-mantle boundary is unknown, it generates a source of uncertainties, and may influence model results.More detail on error sources and model uncertainties/limitations can be found in Ismail-Zadeh & Tackley50 .While the inclusion of the model refinements and reduction of uncertainties are worthwhile, the numerical results derived from this simple model illustrate how the thermal diffusion influences the mantle plume evolution.
and high-resolution geophysical data acquisition, the detailed dynamic reconstruction of a prominent mantle plume from its present diffused state should become possible in the near feature.A source of errors can arise from the geometry of a model domain, for example, a three-dimensional Cartesian model versus a spherical shell model.A detailed comparison of both models of the thermal convection showed that the mean temperature in each model attains almost the same value at the dimensionless domain size of 3×3×1 and at Rayleigh numbers greater than 10 5 (ref.53).The dimensionless size of the model domain used in this study is

Table S1 .
Euler poles describing the motion of Eurasia with respect to the mantle.